# This file replicates the analysis and graphics in
# Wibbels and Ahlquist 2011 ISQ
# Last updated: 1-12-2011						
###################################################
# Important Comments
# ------------------
# Open R by double-clicking on this file to 
#   automatically set the appropriate working
#   directory.
# ------------------
# If the packages below are not already installed,
#   use the syntax install.packages("foreign"), to 
#   for example, install the foreign package.
# ------------------

# Create useful functions

na.mean<-function(x){
    m<-mean(x,na.rm=T)
    return(m)
    }

na.sd<-function(x){
    m<-sd(x,na.rm=T)
    return(m)
    }

na.coef.var.restrict<-function(x){
    if(length(na.omit(x))>2) na.sd(x)/abs(na.mean(x))
    else NA
    }

# Load required packages

library(foreign)
library(MASS)
library(car)
library(systemfit)
library(spdep)
library(lattice)

# Load and subset data

wd<-read.csv("WibbelsAhlquistISQ2011data.csv")

wd2<-subset(wd,(dev==1) & (zaireweird==0))
wd1<-subset(wd,(dev==1) & (excom==0) & (zaireweird==0))

wd.1982<-subset(wd1, year==1982)
wd.1995<-subset(wd1, year==1995)

# Figure 1

data.early<-subset(wd1,(year<1983)&(dev==1)&(excom==0)&(zaireweird==0))
ssex.earl<-cbind(data.early$soc.sec.exp.gdp.final,data.early$soc.sec.exp.govex.final)
row.names(ssex.earl)<-data.early$wbcode
stability<-aggregate(ssex.earl,list(data.early$wbcode),na.coef.var.restrict)
rownames(stability)<-stability$Group.1

boxplot(stability[,2:3],names=c("Soc.Sec/GDP", "Soc.Sec/GovEx"),
    ylab="Coef. of variation", boxwex=.2, staplewex=.1, notch=T)

# Figure 2

palette(c("black", "black"))
scatterplot(wd.1982$soc.sec.exp.govex.final.early ~ I(wd.1982$isi1early+51),
span=.5, xlab="ISI pre-83 avg", ylab="Soc. sec. % gov. spend pre-83 avg", xlim=c(20,85))   

# Set up variables for Model 1

mdd<-na.omit(read.csv("mdd1982.csv"))

outs<-c("NAM", "ABW", "ATG", "BHS", "BMU", "BRB", "COM", "CPV", "DMA", "ERI", "GRD", "JAM", "KNA", "KIR", "LCA", "MAC", "MDG", "MUS", "STP", "SYC", "TON", "VCT", "VUT")
ins<-setdiff(wd.1982$wbcode,outs)
wd.1982.r<-wd.1982[wd.1982$wbcode%in%ins,]
  
nblist <- vector(mode="list",length=dim(wd.1982.r)[1])
attr(nblist,"region.id") <- wd.1982.r$wbcode
attr(nblist,"class") <- "nb"
nbnms <- data.frame(wd.1982.r$wbcode,c(1:dim(wd.1982.r)[1]))
names(nbnms) <-c ("acr","nm")
min100<-mdd

nodata <- setdiff(wd.1982.r$wbcode,unique(c(as.character(min100$wba),as.character(min100$wbb))))
# Find neighbors for each row in the sldv
for(i in 1:dim(wd.1982.r)[1]){
    temp <- min100[min100$wba==as.character(wd.1982.r$wbcode[i]) |
    min100$wbb==as.character(wd.1982.r$wbcode[i]),]
    cty <- unique(c(as.character(temp$wba),as.character(temp$wbb)))
    cty <- setdiff(cty,wd.1982$wbcode.r[i])
    nblist[[i]] <- nbnms[match(cty,nbnms$acr),"nm"]
}

# Model 1 (Table 1)

model.1<-lm(isi1early~marketsizeearly + laborendow_predebt + familyfarmpct_predebt + resourceexppredebt, data=wd.1982)
summary(model.1)
sqrt(diag(hccm(model.1)))

model.1.sl <- lagsarlm(isi1early~marketsizeearly + laborendow_predebt + familyfarmpct_predebt + resourceexppredebt, data=wd.1982.r, nb2listw(nblist, zero.policy=T), method="eigen", quiet=T, na.action = na.omit)
summary(model.1.sl)

# Set up variables for Models 2 and 3

isi.hat.sl<-predict(model.1.sl)
names(isi.hat.sl)<-wd.1982.r[names(isi.hat.sl),"wbcode"]
names(isi.hat.sl)<-match(names(isi.hat.sl), wd.1982$wbcode)
isi.hats.sl<-isi.hats<-rep(NA, nrow(wd.1982))
isi.hats.sl[as.numeric(names(isi.hat.sl))]<-isi.hat.sl

# Models 2 and 3 (Table 2)

model.2<- lm(soc.sec.exp.govex.final.early~ isi1early + popoldearly + politypredebt + log(rgdppcpppchpwt_predebt), data=wd.1982)
summary(model.2)
sqrt(diag(hccm(model.2)))

model.3<- lm(soc.sec.exp.govex.final.early~ isi.hats.sl + popoldearly + politypredebt + log(rgdppcpppchpwt_predebt), data=wd.1982)
summary(model.3)
sqrt(diag(hccm(model.3)))